% setup env
close all;
clear all;
clc;

% Probe Signal Direction
alpha = deg2rad(-180:10:180);
beta = deg2rad(-90:10:90);
%index for b0
b0 = ceil(length(beta)/2);

% Signal direction
alpha_s = deg2rad(20); %30);
beta_s = deg2rad(0);

% Jam Direction
alpha_j = deg2rad(30:5:95);
beta_j = 0;

% Diagram of Single Element 
Fa = ones(length(beta),1) * (1+cos(alpha-pi/2)) .^2;
[alpha_mesh,beta_mesh] = meshgrid(alpha, beta);
[Fax, Fay, Faz] = sph2cart(beta_mesh, alpha_mesh,Fa);

% Plot 
figure;
surf(Fax, Fay, Faz);
title('Diagram of Single Element');
xlabel('X');
ylabel('Y');
zlabel('Z');
axis equal;

% jam/noise
q_db = -10; %10; % in db
q = 10^(q_db / 10); % absolute

% Allocate memory
F = nan(length(beta), length(alpha));

figure(2);
for alpha_j_cur = deg2rad(70) %alpha_j
    % Compute beta_w
    H_j = H(alpha_j_cur, beta_j);
    D = q * H_j * H_j' + eye(2);
    H_s = H(alpha_s, beta_s);
    % A \ B => INV(A)*B
    beta_w = (D \ H_s) / (H_s' * ( D \ H_s ));
    for a = 1:length(alpha)
        for b = 1:length(beta)
            % Compute Diagram
            F(b, a) = abs(beta_w.' * H(alpha(a) ,beta(b))) .^2;
        end
    end
    % Correction Diagram by Single Element
    F = F .* Fa;
    
    % Plotting
    [Fax, Fay, Faz] = sph2cart(beta_mesh, alpha_mesh, F);
    % Plot 
    %figure;
    subplot(1,2,1)
    surf(Fax, Fay, Faz);
    title('Diagram');
    xlabel('X');
    ylabel('Y');
    zlabel('Z');
    axis equal;
    
    % Cut Plane
    subplot(1,2,2)
    Fcp = F(b0,:);
    polar(alpha, Fcp);
    title(sprintf('Jam a=%0.2f, Signal a=%0.2f',rad2deg(alpha_j_cur),rad2deg(alpha_s)));
    hold on
    polar([alpha_s alpha_s], [0 max(Fcp)],'green');
    polar([alpha_j_cur alpha_j_cur], [0 max(Fcp)],'red');
    hold off
    
    filename = sprintf('pic/a_j=%03.0f.png',rad2deg(alpha_j_cur));
    saveas(2,filename);
end

%
